/* -*- mode: C -*-  */
/* 
   IGraph library.
   Copyright (C) 2003-2012  Gabor Csardi <csardi.gabor@gmail.com>
   334 Harvard street, Cambridge, MA 02139 USA
   
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.
   
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
   
   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 
   02110-1301 USA

*/

#include "igraph_memory.h"
#include "igraph_error.h"
#include "igraph_random.h"
#include "igraph_qsort.h"

#include <assert.h>
#include <string.h> 		/* memcpy & co. */
#include <stdlib.h>
#include <stdarg.h>		/* va_start & co */
#include <math.h>

/**
 * \ingroup vector
 * \section about_igraph_vector_t_objects About \type igraph_vector_t objects
 * 
 * <para>The \type igraph_vector_t data type is a simple and efficient
 * interface to arrays containing numbers. It is something
 * similar as (but much simpler than) the \type vector template
 * in the C++ standard library.</para>
 *
 * <para>Vectors are used extensively in \a igraph, all
 * functions which expect or return a list of numbers use
 * igraph_vector_t to achieve this.</para>
 *
 * <para>The \type igraph_vector_t type usually uses
 * O(n) space
 * to store n elements. Sometimes it
 * uses more, this is because vectors can shrink, but even if they
 * shrink, the current implementation does not free a single bit of
 * memory.</para>
 * 
 * <para>The elements in an \type igraph_vector_t
 * object are indexed from zero, we follow the usual C convention
 * here.</para>
 * 
 * <para>The elements of a vector always occupy a single block of
 * memory, the starting address of this memory block can be queried
 * with the \ref VECTOR macro. This way, vector objects can be used
 * with standard mathematical libraries, like the GNU Scientific
 * Library.</para>
 */

/**
 * \ingroup vector
 * \section igraph_vector_constructors_and_destructors Constructors and
 * Destructors
 * 
 * <para>\type igraph_vector_t objects have to be initialized before using
 * them, this is analogous to calling a constructor on them. There are a
 * number of \type igraph_vector_t constructors, for your
 * convenience. \ref igraph_vector_init() is the basic constructor, it
 * creates a vector of the given length, filled with zeros.
 * \ref igraph_vector_copy() creates a new identical copy
 * of an already existing and initialized vector. \ref
 * igraph_vector_init_copy() creates a vector by copying a regular C array. 
 * \ref igraph_vector_init_seq() creates a vector containing a regular
 * sequence with increment one.</para>
 * 
 * <para>\ref igraph_vector_view() is a special constructor, it allows you to
 * handle a regular C array as a \type vector without copying
 * its elements.
 * </para> 
 *
 * <para>If a \type igraph_vector_t object is not needed any more, it
 * should be destroyed to free its allocated memory by calling the
 * \type igraph_vector_t destructor, \ref igraph_vector_destroy().</para>
 * 
 * <para> Note that vectors created by \ref igraph_vector_view() are special,
 * you mustn't call \ref igraph_vector_destroy() on these.</para>
 */

/**
 * \ingroup vector
 * \function igraph_vector_init
 * \brief Initializes a vector object (constructor).
 * 
 * </para><para>
 * Every vector needs to be initialized before it can be used, and
 * there are a number of initialization functions or otherwise called
 * constructors. This function constructs a vector of the given size and
 * initializes each entry to 0. Note that \ref igraph_vector_null() can be
 * used to set each element of a vector to zero. However, if you want a
 * vector of zeros, it is much faster to use this function than to create a
 * vector and then invoke \ref igraph_vector_null().
 * 
 * </para><para>
 * Every vector object initialized by this function should be
 * destroyed (ie. the memory allocated for it should be freed) when it
 * is not needed anymore, the \ref igraph_vector_destroy() function is
 * responsible for this.
 * \param v Pointer to a not yet initialized vector object.
 * \param size The size of the vector.
 * \return error code:
 *       \c IGRAPH_ENOMEM if there is not enough memory.
 * 
 * Time complexity: operating system dependent, the amount of 
 * \quote time \endquote required to allocate
 * O(n) elements,
 * n is the number of elements. 
 */

int FUNCTION(igraph_vector,init)      (TYPE(igraph_vector)* v, int long size) {	
        long int alloc_size= size > 0 ? size : 1;
	if (size < 0) { size=0; }
	v->stor_begin=igraph_Calloc(alloc_size, BASE);
	if (v->stor_begin==0) {
	  IGRAPH_ERROR("cannot init vector", IGRAPH_ENOMEM);
	}
	v->stor_end=v->stor_begin + alloc_size;
	v->end=v->stor_begin+size;

	return 0;
}

/**
 * \ingroup vector
 * \function igraph_vector_view
 * \brief Handle a regular C array as a \type igraph_vector_t.
 * 
 * </para><para>
 * This is a special \type igraph_vector_t constructor. It allows to
 * handle a regular C array as a \type igraph_vector_t temporarily.
 * Be sure that you \em don't ever call the destructor (\ref
 * igraph_vector_destroy()) on objects created by this constructor.
 * \param v Pointer to an uninitialized \type igraph_vector_t object.
 * \param data Pointer, the C array.
 * \param length The length of the C array.
 * \return Pointer to the vector object, the same as the 
 *     \p v parameter, for convenience.
 * 
 * Time complexity: O(1)
 */ 

const TYPE(igraph_vector)*FUNCTION(igraph_vector,view) (const TYPE(igraph_vector) *v,
							const BASE *data, 
							long int length) {
  TYPE(igraph_vector) *v2=(TYPE(igraph_vector)*)v;
  v2->stor_begin=(BASE*)data;
  v2->stor_end=(BASE*)data+length;
  v2->end=v2->stor_end;
  return v;
}

#ifndef BASE_COMPLEX

/**
 * \ingroup vector
 * \function igraph_vector_init_real
 * \brief Create an \type igraph_vector_t from the parameters.
 * 
 * </para><para>
 * Because of how C and the C library handles variable length argument
 * lists, it is required that you supply real constants to this
 * function. This means that
 * \verbatim igraph_vector_t v;
 * igraph_vector_init_real(&amp;v, 5, 1,2,3,4,5); \endverbatim
 * is an error at runtime and the results are undefined. This is
 * the proper way:
 * \verbatim igraph_vector_t v;
 * igraph_vector_init_real(&amp;v, 5, 1.0,2.0,3.0,4.0,5.0); \endverbatim
 * \param v Pointer to an uninitialized \type igraph_vector_t object.
 * \param no Positive integer, the number of \type igraph_real_t
 *    parameters to follow.
 * \param ... The elements of the vector.
 * \return Error code, this can be \c IGRAPH_ENOMEM
 *     if there isn't enough memory to allocate the vector.
 *
 * \sa \ref igraph_vector_init_real_end(), \ref igraph_vector_init_int() for similar
 * functions.
 *
 * Time complexity: depends on the time required to allocate memory,
 * but at least O(n), the number of
 * elements in the vector.
 */

int FUNCTION(igraph_vector,init_real)(TYPE(igraph_vector) *v, int no, ...) {
  int i=0;
  va_list ap;
  IGRAPH_CHECK(FUNCTION(igraph_vector,init)(v, no));

  va_start(ap, no);
  for (i=0; i<no; i++) {
    VECTOR(*v)[i]=(BASE) va_arg(ap, double);
  }
  va_end(ap);
  
  return 0;
}

/**
 * \ingroup vector
 * \function igraph_vector_init_real_end
 * \brief Create an \type igraph_vector_t from the parameters.
 * 
 * </para><para>
 * This constructor is similar to \ref igraph_vector_init_real(), the only
 * difference is that instead of giving the number of elements in the
 * vector, a special marker element follows the last real vector
 * element.
 * \param v Pointer to an uninitialized \type igraph_vector_t object.
 * \param endmark This element will signal the end of the vector. It
 *    will \em not be part of the vector.
 * \param ... The elements of the vector.
 * \return Error code, \c IGRAPH_ENOMEM if there
 *    isn't enough memory.
 * 
 * \sa \ref igraph_vector_init_real() and \ref igraph_vector_init_int_end() for
 * similar functions.
 * 
 * Time complexity: at least O(n) for 
 * n elements plus the time
 * complexity of the memory allocation.
 */

int FUNCTION(igraph_vector,init_real_end)(TYPE(igraph_vector) *v, 
					  BASE endmark, ...) {
  int i=0, n=0;
  va_list ap;

  va_start(ap, endmark);
  while (1) {
    BASE num = (BASE) va_arg(ap, double);
    if (num == endmark) {
      break;
    }
    n++;
  }
  va_end(ap);

  IGRAPH_CHECK(FUNCTION(igraph_vector,init)(v,n));
  IGRAPH_FINALLY(FUNCTION(igraph_vector,destroy), v);
  
  va_start(ap, endmark);
  for (i=0; i<n; i++) {
    VECTOR(*v)[i]=(BASE) va_arg(ap, double);
  }
  va_end(ap);
  
  IGRAPH_FINALLY_CLEAN(1);
  return 0;
}

/**
 * \ingroup vector
 * \function igraph_vector_init_int
 * \brief Create an \type igraph_vector_t containing the parameters.
 * 
 * </para><para>
 * This function is similar to \ref igraph_vector_init_real(), but it expects 
 * \type int parameters. It is important that all parameters
 * should be of this type, otherwise the result of the function call
 * is undefined.
 * \param v Pointer to an uninitialized \type igraph_vector_t object.
 * \param no The number of \type int parameters to follow.
 * \param ... The elements of the vector.
 * \return Error code, \c IGRAPH_ENOMEM if there is
 *    not enough memory.
 * \sa \ref igraph_vector_init_real() and igraph_vector_init_int_end(), these are
 *    similar functions.
 *
 * Time complexity: at least O(n) for 
 * n elements plus the time
 * complexity of the memory allocation.
 */

int FUNCTION(igraph_vector,init_int)(TYPE(igraph_vector) *v, int no, ...) {
  int i=0;
  va_list ap;
  IGRAPH_CHECK(FUNCTION(igraph_vector,init)(v, no));

  va_start(ap, no);
  for (i=0; i<no; i++) {
    VECTOR(*v)[i]=(BASE) va_arg(ap, int);
  }
  va_end(ap);
  
  return 0;
}

/**
 * \ingroup vector
 * \function igraph_vector_init_int_end
 * \brief Create an \type igraph_vector_t from the parameters.
 * 
 * </para><para>
 * This constructor is similar to \ref igraph_vector_init_int(), the only
 * difference is that instead of giving the number of elements in the
 * vector, a special marker element follows the last real vector
 * element.
 * \param v Pointer to an uninitialized \type igraph_vector_t object.
 * \param endmark This element will signal the end of the vector. It
 *    will \em not be part of the vector.
 * \param ... The elements of the vector.
 * \return Error code, \c IGRAPH_ENOMEM if there
 *    isn't enough memory.
 * 
 * \sa \ref igraph_vector_init_int() and \ref igraph_vector_init_real_end() for
 * similar functions.
 *
 * Time complexity: at least O(n) for 
 * n elements plus the time
 * complexity of the memory allocation.
 */

int FUNCTION(igraph_vector_init,int_end)(TYPE(igraph_vector) *v, int endmark, ...) {
  int i=0, n=0;
  va_list ap;

  va_start(ap, endmark);
  while (1) {
    int num = va_arg(ap, int);
    if (num == endmark) {
      break;
    }
    n++;
  }
  va_end(ap);

  IGRAPH_CHECK(FUNCTION(igraph_vector,init)(v, n));
  IGRAPH_FINALLY(FUNCTION(igraph_vector,destroy), v);
  
  va_start(ap, endmark);
  for (i=0; i<n; i++) {
    VECTOR(*v)[i]=(BASE) va_arg(ap, int);
  }
  va_end(ap);
  
  IGRAPH_FINALLY_CLEAN(1);
  return 0;
}

#endif /* ifndef BASE_COMPLEX */

/**
 * \ingroup vector
 * \function igraph_vector_destroy
 * \brief Destroys a vector object.
 *
 * </para><para>
 * All vectors initialized by \ref igraph_vector_init() should be properly
 * destroyed by this function. A destroyed vector needs to be
 * reinitialized by \ref igraph_vector_init(), \ref igraph_vector_init_copy() or
 * another constructor.
 * \param v Pointer to the (previously initialized) vector object to
 *        destroy. 
 *
 * Time complexity: operating system dependent.
 */

void FUNCTION(igraph_vector,destroy)   (TYPE(igraph_vector)* v) {
  assert(v != 0);
  if (v->stor_begin != 0) {
    igraph_Free(v->stor_begin);
    v->stor_begin = NULL;
  }
}

/**
 * \ingroup vector
 * \function igraph_vector_capacity
 * \brief Returns the allocated capacity of the vector
 * 
 * Note that this might be different from the size of the vector (as 
 * queried by \ref igraph_vector_size(), and specifies how many elements 
 * the vector can hold, without reallocation.
 * \param v Pointer to the (previously initialized) vector object
 *          to query.
 * \return The allocated capacity.
 * 
 * \sa \ref igraph_vector_size().
 * 
 * Time complexity: O(1).
 */

long int FUNCTION(igraph_vector,capacity)(const TYPE(igraph_vector)*v) {
  return v->stor_end - v->stor_begin;
}

/**
 * \ingroup vector
 * \function igraph_vector_reserve
 * \brief Reserves memory for a vector.
 * 
 * </para><para>
 * \a igraph vectors are flexible, they can grow and
 * shrink. Growing 
 * however occasionally needs the data in the vector to be copied.
 * In order to avoid this, you can call this function to reserve space for
 * future growth of the vector. 
 * 
 * </para><para>
 * Note that this function does \em not change the size of the
 * vector. Let us see a small example to clarify things: if you
 * reserve space for 100 elements and the size of your
 * vector was (and still is) 60, then you can surely add additional 40
 * elements to your vector before it will be copied.
 * \param v The vector object.
 * \param size The new \em allocated size of the vector.
 * \return Error code:
 *         \c IGRAPH_ENOMEM if there is not enough memory.
 *
 * Time complexity: operating system dependent, should be around
 * O(n), n 
 * is the new allocated size of the vector.
 */

int FUNCTION(igraph_vector,reserve)   (TYPE(igraph_vector)* v, long int size) {
	long int actual_size=FUNCTION(igraph_vector,size)(v);
	BASE *tmp;
	assert(v != NULL);
	assert(v->stor_begin != NULL);
	if (size <= FUNCTION(igraph_vector,size)(v)) { return 0; }

	tmp=igraph_Realloc(v->stor_begin, (size_t) size, BASE);
	if (tmp==0) {
	  IGRAPH_ERROR("cannot reserve space for vector", IGRAPH_ENOMEM);
	}
	v->stor_begin=tmp;
	v->stor_end=v->stor_begin + size;
	v->end=v->stor_begin+actual_size;
	
	return 0;
}

/**
 * \ingroup vector
 * \function igraph_vector_empty
 * \brief Decides whether the size of the vector is zero.
 *
 * \param v The vector object.
 * \return Non-zero number (true) if the size of the vector is zero and
 *         zero (false) otherwise.
 * 
 * Time complexity: O(1).
 */

igraph_bool_t FUNCTION(igraph_vector,empty)     (const TYPE(igraph_vector)* v) {
	assert(v != NULL);
	assert(v->stor_begin != NULL);
	return v->stor_begin == v->end;
}

/**
 * \ingroup vector
 * \function igraph_vector_size
 * \brief Gives the size (=length) of the vector.
 * 
 * \param v The vector object
 * \return The size of the vector.
 *
 * Time complexity: O(1). 
 */

long int FUNCTION(igraph_vector,size)      (const TYPE(igraph_vector)* v) {
	assert(v != NULL);
	assert(v->stor_begin != NULL);
	return v->end - v->stor_begin;
}

/**
 * \ingroup vector
 * \function igraph_vector_clear
 * \brief Removes all elements from a vector.
 * 
 * </para><para>
 * This function simply sets the size of the vector to zero, it does
 * not free any allocated memory. For that you have to call
 * \ref igraph_vector_destroy().
 * \param v The vector object.
 * 
 * Time complexity: O(1).
 */

void FUNCTION(igraph_vector,clear)     (TYPE(igraph_vector)* v) {
	assert(v != NULL);
	assert(v->stor_begin != NULL);
	v->end = v->stor_begin;
}

/**
 * \ingroup vector
 * \function igraph_vector_push_back
 * \brief Appends one element to a vector.
 * 
 * </para><para>
 * This function resizes the vector to be one element longer and
 * sets the very last element in the vector to \p e.
 * \param v The vector object.
 * \param e The element to append to the vector.
 * \return Error code:
 *         \c IGRAPH_ENOMEM: not enough memory.
 * 
 * Time complexity: operating system dependent. What is important is that
 * a sequence of n
 * subsequent calls to this function has time complexity
 * O(n), even if there 
 * hadn't been any space reserved for the new elements by
 * \ref igraph_vector_reserve(). This is implemented by a trick similar to the C++
 * \type vector class: each time more memory is allocated for a
 * vector, the size of the additionally allocated memory is the same
 * as the vector's current length. (We assume here that the time
 * complexity of memory allocation is at most linear.)
 */

int FUNCTION(igraph_vector,push_back) (TYPE(igraph_vector)* v, BASE e) {
  	assert(v != NULL);
	assert(v->stor_begin != NULL);
	
	/* full, allocate more storage */
	if (v->stor_end == v->end) {
		long int new_size = FUNCTION(igraph_vector,size)(v) * 2;
		if (new_size == 0) { new_size = 1; }
		IGRAPH_CHECK(FUNCTION(igraph_vector,reserve)(v, new_size));
	}
	
	*(v->end) = e;
	v->end += 1;
	
	return 0;
}

/**
 * \ingroup vector
 * \function igraph_vector_insert
 * \brief Inserts a single element into a vector.
 *
 * Note that this function does not do range checking. Insertion will shift the
 * elements from the position given to the end of the vector one position to the
 * right, and the new element will be inserted in the empty space created at
 * the given position. The size of the vector will increase by one.
 *
 * \param v The vector object.
 * \param pos The position where the new element is to be inserted.
 * \param value The new element to be inserted.
 */
int FUNCTION(igraph_vector,insert)(TYPE(igraph_vector) *v, long int pos,
				   BASE value) {
  size_t size = (size_t) FUNCTION(igraph_vector,size)(v);
  IGRAPH_CHECK(FUNCTION(igraph_vector,resize)(v, (long) size+1));
  if (pos<size) {
    memmove(v->stor_begin+pos+1, v->stor_begin+pos, 
	    sizeof(BASE)*(size - (size_t) pos));
  }
  v->stor_begin[pos] = value;
  return 0;
}

/**
 * \ingroup vector
 * \section igraph_vector_accessing_elements Accessing elements
 * 
 * <para>The simplest way to access an element of a vector is to use the
 * \ref VECTOR macro. This macro can be used both for querying and setting
 * \type igraph_vector_t elements. If you need a function, \ref
 * igraph_vector_e() queries and \ref igraph_vector_set() sets an element of a
 * vector. \ref igraph_vector_e_ptr() returns the address of an element.</para>
 * 
 * <para>\ref igraph_vector_tail() returns the last element of a non-empty
 * vector. There is no <function>igraph_vector_head()</function> function
 * however, as it is easy to write <code>VECTOR(v)[0]</code>
 * instead.</para>
 */

/**
 * \ingroup vector
 * \function igraph_vector_e
 * \brief Access an element of a vector.
 * \param v The \type igraph_vector_t object.
 * \param pos The position of the element, the index of the first
 *    element is zero.
 * \return The desired element.
 * \sa \ref igraph_vector_e_ptr() and the \ref VECTOR macro.
 * 
 * Time complexity: O(1).
 */

BASE FUNCTION(igraph_vector,e)         (const TYPE(igraph_vector)* v, long int pos) {
	assert(v != NULL);
	assert(v->stor_begin != NULL);
	return * (v->stor_begin + pos);
}

/**
 * \ingroup vector
 * \function igraph_vector_e_ptr
 * \brief Get the address of an element of a vector
 * \param v The \type igraph_vector_t object.
 * \param pos The position of the element, the position of the first
 *   element is zero.
 * \return Pointer to the desired element.
 * \sa \ref igraph_vector_e() and the \ref VECTOR macro.
 * 
 * Time complexity: O(1).
 */

BASE* FUNCTION(igraph_vector,e_ptr)  (const TYPE(igraph_vector)* v, long int pos) {
  assert(v!=NULL);
  assert(v->stor_begin != NULL);
  return v->stor_begin+pos;
}

/**
 * \ingroup vector
 * \function igraph_vector_set
 * \brief Assignment to an element of a vector.
 * \param v The \type igraph_vector_t element.
 * \param pos Position of the element to set.
 * \param value New value of the element.
 * \sa \ref igraph_vector_e().
 */

void FUNCTION(igraph_vector,set)       (TYPE(igraph_vector)* v, 
					long int pos, BASE value) {
	assert(v != NULL);
	assert(v->stor_begin != NULL);	
	*(v->stor_begin + pos) = value;
}

/**
 * \ingroup vector
 * \function igraph_vector_null
 * \brief Sets each element in the vector to zero.
 * 
 * </para><para>
 * Note that \ref igraph_vector_init() sets the elements to zero as well, so
 * it makes no sense to call this function on a just initialized
 * vector. Thus if you want to construct a vector of zeros, then you should
 * use \ref igraph_vector_init().
 * \param v The vector object.
 *
 * Time complexity: O(n), the size of
 * the vector. 
 */

void FUNCTION(igraph_vector,null)      (TYPE(igraph_vector)* v) {
	assert(v != NULL);
	assert(v->stor_begin != NULL);
	if (FUNCTION(igraph_vector,size)(v)>0) {
		memset(v->stor_begin, 0, 
		       sizeof(BASE)*(size_t) FUNCTION(igraph_vector,size)(v));
	}
}

/**
 * \function igraph_vector_fill
 * \brief Fill a vector with a constant element
 * 
 * Sets each element of the vector to the supplied constant.
 * \param vector The vector to work on.
 * \param e The element to fill with.
 * 
 * Time complexity: O(n), the size of the vector.
 */

void FUNCTION(igraph_vector,fill)      (TYPE(igraph_vector)* v, BASE e) {
  BASE *ptr;
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  for (ptr = v->stor_begin; ptr < v->end; ptr++) {
    *ptr = e;
  }  
}  

/**
 * \ingroup vector
 * \function igraph_vector_tail
 * \brief Returns the last element in a vector.
 *
 * </para><para>
 * It is an error to call this function on an empty vector, the result
 * is undefined.
 * \param v The vector object.
 * \return The last element.
 * 
 * Time complexity: O(1).
 */

BASE FUNCTION(igraph_vector,tail)(const TYPE(igraph_vector) *v) {
  assert(v!=NULL);
  assert(v->stor_begin != NULL);
  return *((v->end)-1);
}

/**
 * \ingroup vector
 * \function igraph_vector_pop_back
 * \brief Removes and returns the last element of a vector.
 *
 * </para><para>
 * It is an error to call this function with an empty vector.
 * \param v The vector object.
 * \return The removed last element.
 * 
 * Time complexity: O(1).
 */

BASE FUNCTION(igraph_vector,pop_back)(TYPE(igraph_vector)* v) {
  BASE tmp;
  assert(v!=NULL);
  assert(v->stor_begin != NULL);
  assert(v->end != v->stor_begin);
  tmp=FUNCTION(igraph_vector,e)(v, FUNCTION(igraph_vector,size)(v)-1);
  v->end -= 1;
  return tmp;
}

#ifndef NOTORDERED

/**
 * \ingroup vector
 * \function igraph_vector_sort_cmp
 * \brief Internal comparison function of vector elements, used by 
 * \ref igraph_vector_sort().
 */

int FUNCTION(igraph_vector,sort_cmp)(const void *a, const void *b) {
  const BASE *da = (const BASE *) a;
  const BASE *db = (const BASE *) b;

  return (*da > *db) - (*da < *db);
}

/**
 * \ingroup vector
 * \function igraph_vector_sort
 * \brief Sorts the elements of the vector into ascending order.
 * 
 * </para><para>
 * This function uses the built-in sort function of the C library.
 * \param v Pointer to an initialized vector object.
 *
 * Time complexity: should be
 * O(nlogn) for
 * n 
 * elements.
 */

void FUNCTION(igraph_vector,sort)(TYPE(igraph_vector) *v) {
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  igraph_qsort(v->stor_begin, (size_t) FUNCTION(igraph_vector,size)(v), 
    sizeof(BASE), FUNCTION(igraph_vector,sort_cmp));
}

/**
 * Ascending comparison function passed to qsort from  igraph_vector_qsort_ind
 */
int FUNCTION(igraph_vector,i_qsort_ind_cmp_asc)(const void *p1 , const void *p2) {
  BASE **pa = (BASE **) p1;
  BASE **pb = (BASE **) p2;
  if( **pa < **pb ) return -1;
  if( **pa > **pb) return 1;
  return 0;
}

/**
 * Descending comparison function passed to qsort from  igraph_vector_qsort_ind
 */
int FUNCTION(igraph_vector,i_qsort_ind_cmp_desc)(const void *p1 , const void *p2) {
  BASE **pa = (BASE **) p1;
  BASE **pb = (BASE **) p2;
  if( **pa < **pb ) return 1;
  if( **pa > **pb) return -1;
  return 0;
}

/**
 * \function igraph_vector_qsort_ind
 * \brief Return a permutation of indices that sorts a vector 
 *
 * Takes an unsorted array \c v as input and computes an array of
 * indices inds such that v[ inds[i] ], with i increasing from 0, is
 * an ordered array (either ascending or descending, depending on
 * \v order). The order of indices for identical elements is not
 * defined.
 *
 * \param v the array to be sorted
 * \param inds the output array of indices. this must be initialized,
 *         but will be resized 
 * \param descending whether the output array should be sorted in descending
 *        order.
 * \return Error code.
 *
 * This routine uses the C library qsort routine.
 * Algorithm: 1) create an array of pointers to the elements of v. 2)
 * Pass this array to qsort. 3) after sorting the difference between
 * the pointer value and the first pointer value gives its original
 * position in the array. Use this to set the values of inds. 
 *
 * Some tests show that this routine is faster than
 * igraph_vector_heapsort_ind by about 10 percent 
 * for small vectors to a factor of two for large vectors.
 */

long int FUNCTION(igraph_vector,qsort_ind)(TYPE(igraph_vector) *v, 
					   igraph_vector_t *inds, igraph_bool_t descending) {
  long int i;
  BASE **vind, *first;
  size_t n = (size_t) FUNCTION(igraph_vector,size)(v);
  IGRAPH_CHECK(igraph_vector_resize(inds, (long) n));
  if (n==0) { return 0; }
  vind = igraph_Calloc(n, BASE*);
  if (vind == 0) {
    IGRAPH_ERROR("igraph_vector_qsort_ind failed", IGRAPH_ENOMEM);
  }
  for(i=0; i<n; i++) {
    vind[i] = &VECTOR(*v)[i];
  }
  first = vind[0];
  if (descending)
    igraph_qsort(vind, n, sizeof(BASE**), FUNCTION(igraph_vector,i_qsort_ind_cmp_desc));
  else
    igraph_qsort(vind, n, sizeof(BASE**), FUNCTION(igraph_vector,i_qsort_ind_cmp_asc));
  for(i=0; i<n; i++) {
    VECTOR(*inds)[i] = vind[i]-first;
  }
  igraph_Free(vind);
  return 0;
}

#endif

/**
 * \ingroup vector
 * \function igraph_vector_resize
 * \brief Resize the vector.
 *
 * </para><para>
 * Note that this function does not free any memory, just sets the
 * size of the vector to the given one. It can on the other hand 
 * allocate more memory if the new size is larger than the previous
 * one. In this case the newly appeared elements in the vector are
 * \em not set to zero, they are uninitialized.
 * \param v The vector object
 * \param newsize The new size of the vector.
 * \return Error code, 
 *         \c IGRAPH_ENOMEM if there is not enough
 *         memory. Note that this function \em never returns an error
 *         if the vector is made smaller.
 * \sa \ref igraph_vector_reserve() for allocating memory for future
 * extensions of a vector. \ref igraph_vector_resize_min() for 
 * deallocating the unnneded memory for a vector.
 * 
 * Time complexity: O(1) if the new
 * size is smaller, operating system dependent if it is larger. In the
 * latter case it is usually around
 * O(n),
 * n is the new size of the vector. 
 */

int FUNCTION(igraph_vector,resize)(TYPE(igraph_vector)* v, long int newsize) {
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  IGRAPH_CHECK(FUNCTION(igraph_vector,reserve)(v, newsize));
  v->end = v->stor_begin+newsize;
  return 0;
}

/**
 * \ingroup vector
 * \function igraph_vector_resize_min
 * \brief Deallocate the unused memory of a vector.
 * 
 * </para><para>
 * Note that this function involves additional memory allocation and 
 * may result an out-of-memory error.
 * \param v Pointer to an initialized vector.
 * \return Error code.
 *
 * \sa \ref igraph_vector_resize(), \ref igraph_vector_reserve().
 * 
 * Time complexity: operating system dependent. 
 */

int FUNCTION(igraph_vector,resize_min)(TYPE(igraph_vector)*v) {
  size_t size;
  BASE *tmp;
  if (v->stor_end == v->end) { return 0; }
  
  size = (size_t) (v->end - v->stor_begin);
  tmp=igraph_Realloc(v->stor_begin, size, BASE);
  if (tmp==0) {
    IGRAPH_ERROR("cannot resize vector", IGRAPH_ENOMEM);
  } else {
    v->stor_begin = tmp;
    v->stor_end = v->end = v->stor_begin + size;
  }
  
  return 0;
}

#ifndef NOTORDERED

/**
 * \ingroup vector
 * \function igraph_vector_max
 * \brief Gives the maximum element of the vector.
 *
 * </para><para>
 * If the size of the vector is zero, an arbitrary number is
 * returned.
 * \param v The vector object.
 * \return The maximum element.
 *
 * Time complexity: O(n),
 * n is the size of the vector. 
 */

BASE FUNCTION(igraph_vector,max)(const TYPE(igraph_vector)* v) {
  BASE max;
  BASE *ptr;
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  max=*(v->stor_begin);
  ptr=v->stor_begin+1;
  while (ptr < v->end) {
    if ((*ptr) > max) {
      max=*ptr;
    }
    ptr++;
  }
  return max;
}

/**
 * \ingroup vector
 * \function igraph_vector_which_max
 * \brief Gives the position of the maximum element of the vector.
 *
 * </para><para>
 * If the size of the vector is zero, -1 is 
 * returned.
 * \param v The vector object.
 * \return The position of the first maximum element.
 *
 * Time complexity: O(n),
 * n is the size of the vector. 
 */

long int FUNCTION(igraph_vector,which_max)(const TYPE(igraph_vector)* v) {
  long int which=-1;
  if (!FUNCTION(igraph_vector,empty)(v)) {
    BASE max;
    BASE *ptr;
    long int pos;
    assert(v != NULL);
    assert(v->stor_begin != NULL);
    max=*(v->stor_begin); which=0;
    ptr=v->stor_begin+1; pos=1;
    while (ptr < v->end) {
      if ((*ptr) > max) {
	max=*ptr;
	which=pos;
      }
      ptr++; pos++;
    }
  }
  return which;
}

/**
 * \function igraph_vector_min
 * \brief Smallest element of a vector.
 * 
 * The vector must be non-empty.
 * \param v The input vector.
 * \return The smallest element of \p v.
 * 
 * Time complexity: O(n), the number of elements.
 */

BASE FUNCTION(igraph_vector,min)(const TYPE(igraph_vector)* v) {
  BASE min;
  BASE *ptr;
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  min=*(v->stor_begin);
  ptr=v->stor_begin+1;
  while (ptr < v->end) {
    if ((*ptr) < min) {
      min=*ptr;
    }
    ptr++;
  }
  return min;
}

/**
 * \function igraph_vector_which_min
 * \brief Index of the smallest element.
 * 
 * The vector must be non-empty.
 * If the smallest element is not unique, then the index of the first
 * is returned.
 * \param v The input vector.
 * \return Index of the smallest element.
 * 
 * Time complexity: O(n), the number of elements.
 */

long int FUNCTION(igraph_vector,which_min)(const TYPE(igraph_vector)* v) {
  long int which=-1;
  if (!FUNCTION(igraph_vector,empty)(v)) {
    BASE min;
    BASE *ptr;
    long int pos;
    assert(v != NULL);
    assert(v->stor_begin != NULL);
    min=*(v->stor_begin); which=0;
    ptr=v->stor_begin+1; pos=1;
    while (ptr < v->end) {
      if ((*ptr) < min) {
	min=*ptr;
	which=pos;
      }
      ptr++; pos++;
    }
  }
  return which;
}

#endif

/**
 * \ingroup vector
 * \function igraph_vector_init_copy
 * \brief Initializes a vector from an ordinary C array (constructor).
 * 
 * \param v Pointer to an uninitialized vector object.
 * \param data A regular C array.
 * \param length The length of the C array.
 * \return Error code: 
 *         \c IGRAPH_ENOMEM if there is not enough memory.
 * 
 * Time complexity: operating system specific, usually
 * O(\p length).
 */

int FUNCTION(igraph_vector,init_copy)(TYPE(igraph_vector) *v, 
				      BASE *data, long int length) {
  v->stor_begin=igraph_Calloc(length, BASE);
  if (v->stor_begin==0) {
    IGRAPH_ERROR("cannot init vector from array", IGRAPH_ENOMEM);
  }
  v->stor_end=v->stor_begin+length;
  v->end=v->stor_end;
  memcpy(v->stor_begin, data, (size_t) length * sizeof(BASE));
  
  return 0;
}

/**
 * \ingroup vector
 * \function igraph_vector_copy_to
 * \brief Copies the contents of a vector to a C array.
 * 
 * </para><para>
 * The C array should have sufficient length.
 * \param v The vector object.
 * \param to The C array.
 * 
 * Time complexity: O(n),
 * n is the size of the vector.
 */

void FUNCTION(igraph_vector,copy_to)(const TYPE(igraph_vector) *v, BASE *to) {
				     
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  if (v->end != v->stor_begin) {
    memcpy(to, v->stor_begin, sizeof(BASE) * (size_t) (v->end - v->stor_begin));
  }
}

/**
 * \ingroup vector
 * \function igraph_vector_copy
 * \brief Initializes a vector from another vector object (constructor).
 * 
 * </para><para>
 * The contents of the existing vector object will be copied to
 * the new one.
 * \param to Pointer to a not yet initialized vector object.
 * \param from The original vector object to copy.
 * \return Error code:
 *         \c IGRAPH_ENOMEM if there is not enough memory.
 * 
 * Time complexity: operating system dependent, usually
 * O(n),
 * n is the size of the vector. 
 */

int FUNCTION(igraph_vector,copy)(TYPE(igraph_vector) *to, 
				 const TYPE(igraph_vector) *from) {
  assert(from != NULL);
  assert(from->stor_begin != NULL);
  to->stor_begin=igraph_Calloc(FUNCTION(igraph_vector,size)(from), BASE);
  if (to->stor_begin==0) {
    IGRAPH_ERROR("cannot copy vector", IGRAPH_ENOMEM);
  }
  to->stor_end=to->stor_begin+FUNCTION(igraph_vector,size)(from);
  to->end=to->stor_end;
  memcpy(to->stor_begin, from->stor_begin, 
	 (size_t) FUNCTION(igraph_vector,size)(from) * sizeof(BASE));
  
  return 0;
}

/**
 * \ingroup vector
 * \function igraph_vector_sum
 * \brief Calculates the sum of the elements in the vector.
 *
 * </para><para>
 * For the empty vector 0.0 is returned.
 * \param v The vector object.
 * \return The sum of the elements.
 * 
 * Time complexity: O(n), the size of
 * the vector. 
 */

BASE FUNCTION(igraph_vector,sum)(const TYPE(igraph_vector) *v) {
  BASE res=ZERO;
  BASE *p;
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  for (p=v->stor_begin; p<v->end; p++) {
#ifdef SUM
    SUM(res,res,*p);
#else
    res += *p;
#endif
  }
  return res;
}

igraph_real_t FUNCTION(igraph_vector,sumsq)(const TYPE(igraph_vector) *v) {
  igraph_real_t res=0.0;
  BASE *p;
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  for (p=v->stor_begin; p<v->end; p++) {
#ifdef SQ
    res += SQ(*p);
#else
    res += (*p) * (*p);
#endif
  }
  return res;
}
  
/**
 * \ingroup vector
 * \function igraph_vector_prod
 * \brief Calculates the product of the elements in the vector.
 * 
 * </para><para>
 * For the empty vector one (1) is returned.
 * \param v The vector object.
 * \return The product of the elements.
 * 
 * Time complexity: O(n), the size of
 * the vector. 
 */

BASE FUNCTION(igraph_vector,prod)(const TYPE(igraph_vector) *v) {
  BASE res=ONE;
  BASE *p;
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  for (p=v->stor_begin; p<v->end; p++) {
#ifdef PROD
    PROD(res,res,*p);
#else
    res *= *p;
#endif
  }
  return res;
}

/**
 * \ingroup vector
 * \function igraph_vector_cumsum
 * \brief Calculates the cumulative sum of the elements in the vector.
 *
 * </para><para>
 * \param to An initialized vector object that will store the cumulative
 *           sums. Element i of this vector will store the sum of the elements
 *           of the 'from' vector, up to and including element i.
 * \param from The input vector.
 * \return Error code.
 * 
 * Time complexity: O(n), the size of the vector. 
 */

int FUNCTION(igraph_vector,cumsum)(TYPE(igraph_vector) *to,
                                   const TYPE(igraph_vector) *from) {
  BASE res=ZERO;
  BASE *p, *p2;

  assert(from != NULL);
  assert(from->stor_begin != NULL);
  assert(to != NULL);
  assert(to->stor_begin != NULL);

  IGRAPH_CHECK(FUNCTION(igraph_vector,resize)(to, FUNCTION(igraph_vector,size)(from)));

  for (p = from->stor_begin, p2 = to->stor_begin; p < from->end; p++, p2++) {
#ifdef SUM
    SUM(res,res,*p);
#else
    res += *p;
#endif
    *p2 = res;
  }

  return 0;
}

#ifndef NOTORDERED

/**
 * \ingroup vector
 * \function igraph_vector_init_seq
 * \brief Initializes a vector with a sequence.
 * 
 * </para><para>
 * The vector will contain the numbers \p from,
 * \p from+1, ..., \p to.
 * \param v Pointer to an uninitialized vector object.
 * \param from The lower limit in the sequence (inclusive).
 * \param to The upper limit in the sequence (inclusive).
 * \return Error code:
 *         \c IGRAPH_ENOMEM: out of memory.
 *
 * Time complexity: O(n), the number
 * of elements in the vector. 
 */

int FUNCTION(igraph_vector,init_seq)(TYPE(igraph_vector) *v, 
				     BASE from, BASE to) {
  BASE *p;
  IGRAPH_CHECK(FUNCTION(igraph_vector,init)(v, (long int) (to-from+1)));

  for (p=v->stor_begin; p<v->end; p++) {
    *p = from++;
  }
  
  return 0;
}

#endif

/**
 * \ingroup vector
 * \function igraph_vector_remove_section
 * \brief Deletes a section from a vector.
 * 
 * </para><para>
 * Note that this function does not do range checking. The result is
 * undefined if you supply invalid limits.
 * \param v The vector object.
 * \param from The position of the first element to remove.
 * \param to The position of the first element \em not to remove.
 *
 * Time complexity: O(n-from),
 * n is the number of elements in the
 * vector. 
 */

void FUNCTION(igraph_vector,remove_section)(TYPE(igraph_vector) *v, 
					    long int from, long int to) {
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  /* Not removing from the end? */
  if (to < FUNCTION(igraph_vector,size)(v)) {
    memmove(v->stor_begin+from, v->stor_begin+to,
	    sizeof(BASE) * (size_t) (v->end-v->stor_begin-to));
  }
  v->end -= (to-from);
}

/**
 * \ingroup vector
 * \function igraph_vector_remove
 * \brief Removes a single element from a vector.
 *
 * Note that this function does not do range checking.
 * \param v The vector object.
 * \param elem The position of the element to remove.
 * 
 * Time complexity: O(n-elem),
 * n is the number of elements in the
 * vector. 
 */

void FUNCTION(igraph_vector,remove)(TYPE(igraph_vector) *v, long int elem) {
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  FUNCTION(igraph_vector,remove_section)(v, elem, elem+1);
}

/**
 * \ingroup vector
 * \function igraph_vector_move_interval
 * \brief Copies a section of a vector.
 *
 * </para><para>
 * The result of this function is undefined if the source and target
 * intervals overlap.
 * \param v The vector object.
 * \param begin The position of the first element to move.
 * \param end The position of the first element \em not to move.
 * \param to The target position.
 * \return Error code, the current implementation always returns with
 *    success. 
 *
 * Time complexity: O(end-begin).
 */

int FUNCTION(igraph_vector,move_interval)(TYPE(igraph_vector) *v, 
					  long int begin, long int end, 
					  long int to) {
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  memcpy(v->stor_begin+to, v->stor_begin+begin, 
	 sizeof(BASE) * (size_t) (end-begin));

  return 0;
}

int FUNCTION(igraph_vector,move_interval2)(TYPE(igraph_vector) *v, 
					  long int begin, long int end, 
					  long int to) {
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  memmove(v->stor_begin+to, v->stor_begin+begin, 
	  sizeof(BASE) * (size_t) (end-begin));

  return 0;
}

/**
 * \ingroup vector
 * \function igraph_vector_permdelete
 * \brief Remove elements of a vector (for internal use).
 */

void FUNCTION(igraph_vector,permdelete)(TYPE(igraph_vector) *v, 
					const igraph_vector_t *index, long int nremove) {
  long int i, n;
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  n = FUNCTION(igraph_vector,size)(v);
  for (i=0; i<n; i++) {
    if (VECTOR(*index)[i] != 0) {
      VECTOR(*v)[ (long int)VECTOR(*index)[i]-1 ] = VECTOR(*v)[i];
    }
  }
  v->end -= nremove;
}

#ifndef NOTORDERED

/**
 * \ingroup vector
 * \function igraph_vector_isininterval
 * \brief Checks if all elements of a vector are in the given
 * interval.
 * 
 * \param v The vector object.
 * \param low The lower limit of the interval (inclusive).
 * \param high The higher limit of the interval (inclusive).
 * \return True (positive integer) if all vector elements are in the
 *   interval, false (zero) otherwise.
 *
 * Time complexity: O(n), the number
 * of elements in the vector.
 */

igraph_bool_t FUNCTION(igraph_vector,isininterval)(const TYPE(igraph_vector) *v, 
						   BASE low, 
						   BASE high) {
  BASE *ptr;
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  for (ptr=v->stor_begin; ptr<v->end; ptr++) {
    if (*ptr < low || *ptr >high) {
      return 0;
    }
  }
  return 1;
}

/**
 * \ingroup vector
 * \function igraph_vector_any_smaller
 * \brief Checks if any element of a vector is smaller than a limit.
 * 
 * \param v The \type igraph_vector_t object.
 * \param limit The limit.
 * \return True (positive integer) if the vector contains at least one
 *   smaller element than \p limit, false (zero)
 *   otherwise. 
 * 
 * Time complexity: O(n), the number
 * of elements in the vector.
 */

igraph_bool_t FUNCTION(igraph_vector,any_smaller)(const TYPE(igraph_vector) *v, 
						  BASE limit) {
  BASE *ptr;
  assert(v != NULL);
  assert(v->stor_begin != NULL);
  for (ptr=v->stor_begin; ptr<v->end; ptr++) {
    if (*ptr < limit) {
      return 1;
    }
  }
  return 0;
}

#endif

/**
 * \ingroup vector
 * \function igraph_vector_all_e
 * \brief Are all elements equal?
 * 
 * \param lhs The first vector.
 * \param rhs The second vector.
 * \return Positive integer (=true) if the elements in the \p lhs are all 
 *    equal to the corresponding elements in \p rhs. Returns \c 0
 *    (=false) if the lengths of the vectors don't match.
 * 
 * Time complexity: O(n), the length of the vectors.
 */

igraph_bool_t FUNCTION(igraph_vector,all_e)(const TYPE(igraph_vector) *lhs, 
					    const TYPE(igraph_vector) *rhs) {
  long int i, s;
  assert(lhs != 0);
  assert(rhs != 0);
  assert(lhs->stor_begin != 0);
  assert(rhs->stor_begin != 0);
  
  s=FUNCTION(igraph_vector,size)(lhs);
  if (s != FUNCTION(igraph_vector,size)(rhs)) {
    return 0;
  } else {
    for (i=0; i<s; i++) {
      BASE l=VECTOR(*lhs)[i];
      BASE r=VECTOR(*rhs)[i];
#ifdef EQ
      if (!EQ(l,r)) {
#else
      if (l!=r) {
#endif
	return 0;
      }
    }
    return 1;
  }
}

igraph_bool_t 
FUNCTION(igraph_vector,is_equal)(const TYPE(igraph_vector) *lhs, 
				 const TYPE(igraph_vector) *rhs) {
  return FUNCTION(igraph_vector,all_e)(lhs, rhs);
}

#ifndef NOTORDERED

/**
 * \ingroup vector
 * \function igraph_vector_all_l
 * \brief Are all elements less?
 * 
 * \param lhs The first vector.
 * \param rhs The second vector.
 * \return Positive integer (=true) if the elements in the \p lhs are all 
 *    less than the corresponding elements in \p rhs. Returns \c 0
 *    (=false) if the lengths of the vectors don't match.
 * 
 * Time complexity: O(n), the length of the vectors.
 */

igraph_bool_t FUNCTION(igraph_vector,all_l)(const TYPE(igraph_vector) *lhs, 
					    const TYPE(igraph_vector) *rhs) {
  long int i, s;
  assert(lhs != 0);
  assert(rhs != 0);
  assert(lhs->stor_begin != 0);
  assert(rhs->stor_begin != 0);
  
  s=FUNCTION(igraph_vector,size)(lhs);
  if (s != FUNCTION(igraph_vector,size)(rhs)) {
    return 0;
  } else {
    for (i=0; i<s; i++) {
      BASE l=VECTOR(*lhs)[i];
      BASE r=VECTOR(*rhs)[i];
      if (l>=r) {
	return 0;
      }
    }
    return 1;
  }  
}

/**
 * \ingroup vector
 * \function igraph_vector_all_g
 * \brief Are all elements greater?
 * 
 * \param lhs The first vector.
 * \param rhs The second vector.
 * \return Positive integer (=true) if the elements in the \p lhs are all 
 *    greater than the corresponding elements in \p rhs. Returns \c 0
 *    (=false) if the lengths of the vectors don't match.
 * 
 * Time complexity: O(n), the length of the vectors.
 */

igraph_bool_t FUNCTION(igraph_vector,all_g)(const TYPE(igraph_vector) *lhs, 
					    const TYPE(igraph_vector) *rhs) {

  long int i, s;
  assert(lhs != 0);
  assert(rhs != 0);
  assert(lhs->stor_begin != 0);
  assert(rhs->stor_begin != 0);
  
  s=FUNCTION(igraph_vector,size)(lhs);
  if (s != FUNCTION(igraph_vector,size)(rhs)) {
    return 0;
  } else {
    for (i=0; i<s; i++) {
      BASE l=VECTOR(*lhs)[i];
      BASE r=VECTOR(*rhs)[i];
      if (l<=r) {
	return 0;
      }
    }
    return 1;
  }  
}

/**
 * \ingroup vector
 * \function igraph_vector_all_le
 * \brief Are all elements less or equal?
 * 
 * \param lhs The first vector.
 * \param rhs The second vector.
 * \return Positive integer (=true) if the elements in the \p lhs are all 
 *    less than or equal to the corresponding elements in \p
 *    rhs. Returns \c 0 (=false) if the lengths of the vectors don't
 *    match.
 * 
 * Time complexity: O(n), the length of the vectors.
 */

igraph_bool_t 
FUNCTION(igraph_vector,all_le)(const TYPE(igraph_vector) *lhs, 
			       const TYPE(igraph_vector) *rhs) {
  long int i, s;
  assert(lhs != 0);
  assert(rhs != 0);
  assert(lhs->stor_begin != 0);
  assert(rhs->stor_begin != 0);
  
  s=FUNCTION(igraph_vector,size)(lhs);
  if (s != FUNCTION(igraph_vector,size)(rhs)) {
    return 0;
  } else {
    for (i=0; i<s; i++) {
      BASE l=VECTOR(*lhs)[i];
      BASE r=VECTOR(*rhs)[i];
      if (l>r) {
	return 0;
      }
    }
    return 1;
  }  
}

/**
 * \ingroup vector
 * \function igraph_vector_all_ge
 * \brief Are all elements greater or equal?
 * 
 * \param lhs The first vector.
 * \param rhs The second vector.
 * \return Positive integer (=true) if the elements in the \p lhs are all 
 *    greater than or equal to the corresponding elements in \p
 *    rhs. Returns \c 0 (=false) if the lengths of the vectors don't
 *    match.
 * 
 * Time complexity: O(n), the length of the vectors.
 */

igraph_bool_t 
FUNCTION(igraph_vector,all_ge)(const TYPE(igraph_vector) *lhs, 
			       const TYPE(igraph_vector) *rhs) {
  long int i, s;
  assert(lhs != 0);
  assert(rhs != 0);
  assert(lhs->stor_begin != 0);
  assert(rhs->stor_begin != 0);
  
  s=FUNCTION(igraph_vector,size)(lhs);
  if (s != FUNCTION(igraph_vector,size)(rhs)) {
    return 0;
  } else {
    for (i=0; i<s; i++) {
      BASE l=VECTOR(*lhs)[i];
      BASE r=VECTOR(*rhs)[i];
      if (l<r) {
	return 0;
      }
    }
    return 1;
  }  
}

#endif

igraph_bool_t FUNCTION(igraph_i_vector,binsearch_slice)(const TYPE(igraph_vector) *v,
                                                        BASE what, long int *pos,
                                                        long int start, long int end);

#ifndef NOTORDERED

/**
 * \ingroup vector
 * \function igraph_vector_binsearch
 * \brief Finds an element by binary searching a sorted vector.
 * 
 * </para><para>
 * It is assumed that the vector is sorted. If the specified element
 * (\p what) is not in the vector, then the
 * position of where it should be inserted (to keep the vector sorted)
 * is returned.
 * \param v The \type igraph_vector_t object.
 * \param what The element to search for.
 * \param pos Pointer to a \type long int. This is set to the
 *   position of an instance of \p what in the
 *   vector if it is present. If \p v does not
 *   contain \p what then
 *   \p pos is set to the position to which it
 *   should be inserted (to keep the the vector sorted of course).
 * \return Positive integer (true) if \p what is
 *   found in the vector, zero (false) otherwise.
 * 
 * Time complexity: O(log(n)),
 * n is the number of elements in
 * \p v.
 */

igraph_bool_t FUNCTION(igraph_vector,binsearch)(const TYPE(igraph_vector) *v, 
						BASE what, long int *pos) {
  return FUNCTION(igraph_i_vector,binsearch_slice)(v, what, pos,
      0, FUNCTION(igraph_vector,size)(v));
}

igraph_bool_t FUNCTION(igraph_i_vector,binsearch_slice)(const TYPE(igraph_vector) *v,
                                                        BASE what, long int *pos,
                                                        long int start, long int end) {
  long int left  = start;
  long int right = end-1;

  while (left <= right) {
    /* (right + left) / 2 could theoretically overflow for long vectors */
    long int middle = left + ((right - left) >> 1);
    if (VECTOR(*v)[middle] > what) {
      right = middle - 1;
    } else if (VECTOR(*v)[middle] < what) {
      left = middle + 1;
    } else {
      if (pos != 0) {
        *pos = middle;
      }
      return 1;
    }
  }

  /* if we are here, the element was not found */
  if (pos != 0) {
    *pos = left;
  }

  return 0;
}

/**
 * \ingroup vector
 * \function igraph_vector_binsearch2
 * \brief Binary search, without returning the index.
 * 
 * </para><para>
 * It is assumed that the vector is sorted.
 * \param v The \type igraph_vector_t object.
 * \param what The element to search for.
 * \return Positive integer (true) if \p what is
 *   found in the vector, zero (false) otherwise.
 * 
 * Time complexity: O(log(n)),
 * n is the number of elements in
 * \p v.
 */

igraph_bool_t FUNCTION(igraph_vector,binsearch2)(const TYPE(igraph_vector) *v, 
						 BASE what) {
  long int left=0;
  long int right=FUNCTION(igraph_vector,size)(v)-1;

  while (left <= right) {
    /* (right + left) / 2 could theoretically overflow for long vectors */
    long int middle = left + ((right - left) >> 1);
    if (what < VECTOR(*v)[middle]) {
      right = middle - 1;
    } else if (what > VECTOR(*v)[middle]) {
      left = middle + 1;
    } else {
      return 1;
    }
  }

  return 0;
}

#endif

/**
 * \function igraph_vector_scale
 * \brief Multiply all elements of a vector by a constant
 * 
 * \param v The vector.
 * \param by The constant.
 * \return Error code. The current implementation always returns with success.
 * 
 * Added in version 0.2.</para><para>
 * 
 * Time complexity: O(n), the number of elements in a vector.
 */

void FUNCTION(igraph_vector,scale)(TYPE(igraph_vector) *v, BASE by) {
  long int i;
  for (i=0; i<FUNCTION(igraph_vector,size)(v); i++) {
#ifdef PROD
    PROD(VECTOR(*v)[i], VECTOR(*v)[i], by);
#else
    VECTOR(*v)[i] *= by;
#endif
  }
}

/**
 * \function igraph_vector_add_constant
 * \brief Add a constant to the vector.
 * 
 * \p plus is added to every element of \p v. Note that overflow
 * might happen.
 * \param v The input vector.
 * \param plus The constant to add.
 * 
 * Time complexity: O(n), the number of elements.
 */

void FUNCTION(igraph_vector,add_constant)(TYPE(igraph_vector) *v, BASE plus) {
  long int i, n=FUNCTION(igraph_vector,size)(v);
  for (i=0; i<n; i++) {
#ifdef SUM
    SUM(VECTOR(*v)[i], VECTOR(*v)[i], plus);
#else
    VECTOR(*v)[i] += plus;
#endif
  }
}

/** 
 * \function igraph_vector_contains
 * \brief Linear search in a vector.
 * 
 * Check whether the supplied element is included in the vector, by
 * linear search.
 * \param v The input vector.
 * \param e The element to look for.
 * \return \c TRUE if the element is found and \c FALSE otherwise.
 * 
 * Time complexity: O(n), the length of the vector.
 */

igraph_bool_t FUNCTION(igraph_vector,contains)(const TYPE(igraph_vector) *v, 
					       BASE e) {
  BASE *p=v->stor_begin;
  while (p<v->end) {
#ifdef EQ
    if (EQ(*p,e)) {
#else
    if (*p==e) { 
#endif
      return 1;
    }
    p++;
  }
  return 0;
}

/**
 * \function igraph_vector_search
 * \brief Search from a given position
 *
 * The supplied element \p what is searched in vector \p v, starting
 * from element index \p from. If found then the index of the first
 * instance (after \p from) is stored in \p pos.
 * \param v The input vector.
 * \param from The index to start searching from. No range checking is
 *     performed. 
 * \param what The element to find.
 * \param pos If not \c NULL then the index of the found element is
 *    stored here.
 * \return Boolean, \c TRUE if the element was found, \c FALSE
 *   otherwise.
 * 
 * Time complexity: O(m), the number of elements to search, the length
 * of the vector minus the \p from argument.
 */

igraph_bool_t FUNCTION(igraph_vector,search)(const TYPE(igraph_vector) *v, 
					     long int from, BASE what, 
					     long int *pos) {
  long int i, n=FUNCTION(igraph_vector,size)(v);  
  for (i=from; i<n; i++) {
#ifdef EQ
    if (EQ(VECTOR(*v)[i], what)) break;
#else
    if (VECTOR(*v)[i]==what) break;
#endif
  }
  
  if (i<n) {
    if (pos != 0) {
      *pos=i;
    }
    return 1;
  } else {
    return 0;
  }
}

#ifndef NOTORDERED

/**
 * \function igraph_vector_filter_smaller
 * \ingroup internal
 */

int FUNCTION(igraph_vector,filter_smaller)(TYPE(igraph_vector) *v, 
					   BASE elem) {
  long int i=0, n=FUNCTION(igraph_vector,size)(v);
  long int s;
  while (i<n && VECTOR(*v)[i]<elem) {
    i++;
  }
  s=i;
  
  while (s<n && VECTOR(*v)[s]==elem) {
    s++;
  }
  
  FUNCTION(igraph_vector,remove_section)(v, 0, i+(s-i)/2);
  return 0;
}

#endif

/**
 * \function igraph_vector_append
 * \brief Append a vector to another one.
 * 
 * The target vector will be resized (except \p from is empty).
 * \param to The vector to append to.
 * \param from The vector to append, it is kept unchanged.
 * \return Error code.
 *
 * Time complexity: O(n), the number of elements in the new vector.
 */

int FUNCTION(igraph_vector,append)(TYPE(igraph_vector) *to, 
				   const TYPE(igraph_vector) *from) {
  long tosize, fromsize;
  
  tosize=FUNCTION(igraph_vector,size)(to);
  fromsize=FUNCTION(igraph_vector,size)(from);
  IGRAPH_CHECK(FUNCTION(igraph_vector,resize)(to, tosize+fromsize));
  memcpy(to->stor_begin+tosize, from->stor_begin, 
	 sizeof(BASE) * (size_t) fromsize);
  to->end=to->stor_begin+tosize+fromsize;
  
  return 0;
}

/**
 * \function igraph_vector_get_interval
 */

int FUNCTION(igraph_vector,get_interval)(const TYPE(igraph_vector) *v, 
					 TYPE(igraph_vector) *res,
					 long int from, long int to) {
  IGRAPH_CHECK(FUNCTION(igraph_vector,resize)(res, to-from));
  memcpy(res->stor_begin, v->stor_begin+from, 
	 (size_t) (to-from) * sizeof(BASE));
  return 0;
}

#ifndef NOTORDERED

/**
 * \function igraph_vector_maxdifference
 * \brief The maximum absolute difference of \p m1 and \p m2
 *
 * The element with the largest absolute value in \p m1 - \p m2 is
 * returned. Both vectors must be non-empty, but they not need to have
 * the same length, the extra elements in the longer vector are ignored.
 * \param m1 The first vector.
 * \param m2 The second vector.
 * \return The maximum absolute difference of \p m1 and \p m2.
 *
 * Time complexity: O(n), the number of elements in the shorter
 * vector.
 */

BASE FUNCTION(igraph_vector,maxdifference)(const TYPE(igraph_vector) *m1,
					   const TYPE(igraph_vector) *m2) {
  long int n1=FUNCTION(igraph_vector,size)(m1);
  long int n2=FUNCTION(igraph_vector,size)(m2);
  long int n= n1 < n2 ? n1 : n2;
  long int i;
  BASE diff=ZERO;
  
  for (i=0; i<n; i++) {
    BASE d=(BASE) fabs(VECTOR(*m1)[i]-VECTOR(*m2)[i]);
    if (d > diff) {
      diff=d;
    }
  }
  
  return diff;
}

#endif

/**
 * \function igraph_vector_update
 * \brief Update a vector from another one.
 * 
 * After this operation the contents of \p to will be exactly the same
 * \p from. \p to will be resized if it was originally shorter or
 * longer than \p from.
 * \param to The vector to update.
 * \param from The vector to update from.
 * \return Error code.
 * 
 * Time complexity: O(n), the number of elements in \p from.
 */

int FUNCTION(igraph_vector,update)(TYPE(igraph_vector) *to, 
				   const TYPE(igraph_vector) *from) {
  size_t n=(size_t) FUNCTION(igraph_vector,size)(from);
  FUNCTION(igraph_vector,resize)(to, (long) n);
  memcpy(to->stor_begin, from->stor_begin, sizeof(BASE)*n);
  return 0;
}

/**
 * \function igraph_vector_swap
 * \brief Swap elements of two vectors.
 * 
 * The two vectors must have the same length, otherwise an error
 * happens.
 * \param v1 The first vector.
 * \param v2 The second vector.
 * \return Error code.
 * 
 * Time complexity: O(n), the length of the vectors.
 */

int FUNCTION(igraph_vector,swap)(TYPE(igraph_vector) *v1, TYPE(igraph_vector) *v2) {
  
  long int i, n1=FUNCTION(igraph_vector,size)(v1);
  long int n2=FUNCTION(igraph_vector,size)(v2);
  if (n1 != n2) {
    IGRAPH_ERROR("Vectors must have the same number of elements for swapping",
		 IGRAPH_EINVAL);
  }
  
  for (i=0; i<n1; i++) {
    BASE tmp;
    tmp=VECTOR(*v1)[i];
    VECTOR(*v1)[i]=VECTOR(*v2)[i];
    VECTOR(*v2)[i]=tmp;
  }
  return 0;
}

/**
 * \function igraph_vector_swap_elements
 * \brief Swap two elements in a vector.
 * 
 * Note that currently no range checking is performed.
 * \param v The input vector.
 * \param i Index of the first element.
 * \param j index of the second element. (Might be the same as the
 * first.)
 * \return Error code, currently always \c IGRAPH_SUCCESS.
 * 
 * Time complexity: O(1).
 */

int FUNCTION(igraph_vector,swap_elements)(TYPE(igraph_vector) *v,
					  long int i, long int j) {
  BASE tmp=VECTOR(*v)[i];
  VECTOR(*v)[i]=VECTOR(*v)[j];
  VECTOR(*v)[j]=tmp;
  
  return 0;
}

/**
 * \function igraph_vector_reverse
 * \brief Reverse the elements of a vector.
 * 
 * The first element will be last, the last element will be
 * first, etc.
 * \param v The input vector.
 * \return Error code, currently always \c IGRAPH_SUCCESS.
 * 
 * Time complexity: O(n), the number of elements.
 */

int FUNCTION(igraph_vector,reverse)(TYPE(igraph_vector) *v) {
  
  long int n=FUNCTION(igraph_vector,size)(v), n2=n/2;
  long int i,j;
  for (i=0, j=n-1; i<n2; i++, j--) {
    BASE tmp;
    tmp=VECTOR(*v)[i];
    VECTOR(*v)[i]=VECTOR(*v)[j];
    VECTOR(*v)[j]=tmp;
  }
  return 0;
}

/**
 * \ingroup vector
 * \function igraph_vector_shuffle
 * \brief Shuffles a vector in-place using the Fisher-Yates method
 * 
 * </para><para>
 * The Fisher-Yates shuffle ensures that every implementation is
 * equally probable when using a proper randomness source. Of course
 * this does not apply to pseudo-random generators as the cycle of
 * these generators is less than the number of possible permutations
 * of the vector if the vector is long enough.
 * \param v The vector object.
 * \return Error code, currently always \c IGRAPH_SUCCESS.
 *
 * Time complexity: O(n),
 * n is the number of elements in the
 * vector. 
 *
 * </para><para>
 * References:
 * \clist
 * \cli (Fisher &amp; Yates 1963)
 *   R. A. Fisher and F. Yates. \emb Statistical Tables for Biological,
 *   Agricultural and Medical Research. \eme Oliver and Boyd, 6th edition,
 *   1963, page 37.
 * \cli (Knuth 1998)
 *   D. E. Knuth. \emb Seminumerical Algorithms, \eme volume 2 of \emb The Art
 *   of Computer Programming. \eme Addison-Wesley, 3rd edition, 1998, page 145.
 * \endclist
 *
 * \example examples/simple/igraph_fisher_yates_shuffle.c
 */

int FUNCTION(igraph_vector,shuffle)(TYPE(igraph_vector) *v) {
  long int n = FUNCTION(igraph_vector,size)(v);
  long int k;
  BASE dummy;

  RNG_BEGIN();
  while (n > 1) {
    k = RNG_INTEGER(0, n-1);
    n--;
    dummy = VECTOR(*v)[n];
    VECTOR(*v)[n] = VECTOR(*v)[k];
    VECTOR(*v)[k] = dummy;
  }
  RNG_END();

  return IGRAPH_SUCCESS;
}

/**
 * \function igraph_vector_add
 * \brief Add two vectors.
 * 
 * Add the elements of \p v2 to \p v1, the result is stored in \p
 * v1. The two vectors must have the same length.
 * \param v1 The first vector, the result will be stored here.
 * \param v2 The second vector, its contents will be unchanged. 
 * \return Error code.
 * 
 * Time complexity: O(n), the number of elements.
 */

int FUNCTION(igraph_vector,add)(TYPE(igraph_vector) *v1, 
				const TYPE(igraph_vector) *v2) {
  
  long int n1=FUNCTION(igraph_vector,size)(v1);
  long int n2=FUNCTION(igraph_vector,size)(v2);
  long int i;
  if (n1 != n2) {
    IGRAPH_ERROR("Vectors must have the same number of elements for swapping",
		 IGRAPH_EINVAL);
  }
  
  for (i=0; i<n1; i++) {
#ifdef SUM
    SUM(VECTOR(*v1)[i], VECTOR(*v1)[i], VECTOR(*v2)[i]);
#else
    VECTOR(*v1)[i] += VECTOR(*v2)[i];
#endif
  }
  
  return 0;
}

/**
 * \function igraph_vector_sub
 * \brief Subtract a vector from another one.
 * 
 * Subtract the elements of \p v2 from \p v1, the result is stored in
 * \p v1. The two vectors must have the same length.
 * \param v1 The first vector, to subtract from. The result is stored
 *    here.
 * \param v2 The vector to subtract, it will be unchanged.
 * \return Error code.
 * 
 * Time complexity: O(n), the length of the vectors.
 */

int FUNCTION(igraph_vector,sub)(TYPE(igraph_vector) *v1, 
				const TYPE(igraph_vector) *v2) {
  
  long int n1=FUNCTION(igraph_vector,size)(v1);
  long int n2=FUNCTION(igraph_vector,size)(v2);
  long int i;
  if (n1 != n2) {
    IGRAPH_ERROR("Vectors must have the same number of elements for swapping",
		 IGRAPH_EINVAL);
  }
  
  for (i=0; i<n1; i++) {
#ifdef DIFF
    DIFF(VECTOR(*v1)[i], VECTOR(*v1)[i], VECTOR(*v2)[i]);
#else
    VECTOR(*v1)[i] -= VECTOR(*v2)[i];
#endif
  }
  
  return 0;
}

/**
 * \function igraph_vector_mul
 * \brief Multiply two vectors.
 * 
 * \p v1 will be multiplied by \p v2, elementwise. The two vectors
 * must have the same length.
 * \param v1 The first vector, the result will be stored here.
 * \param v2 The second vector, it is left unchanged.
 * \return Error code.
 * 
 * Time complexity: O(n), the number of elements.
 */

int FUNCTION(igraph_vector,mul)(TYPE(igraph_vector) *v1, 
				const TYPE(igraph_vector) *v2) {
  
  long int n1=FUNCTION(igraph_vector,size)(v1);
  long int n2=FUNCTION(igraph_vector,size)(v2);
  long int i;
  if (n1 != n2) {
    IGRAPH_ERROR("Vectors must have the same number of elements for swapping",
		 IGRAPH_EINVAL);
  }
  
  for (i=0; i<n1; i++) {
#ifdef PROD
    PROD(VECTOR(*v1)[i], VECTOR(*v1)[i], VECTOR(*v2)[i]);
#else
    VECTOR(*v1)[i] *= VECTOR(*v2)[i];
#endif
  }
  
  return 0;
}

/**
 * \function igraph_vector_div
 * \brief Divide a vector by another one.
 * 
 * \p v1 is divided by \p v2, elementwise. They must have the same length. If the
 * base type of the vector can generate divide by zero errors then
 * please make sure that \p v2 contains no zero if you want to avoid
 * trouble.
 * \param v1 The dividend. The result is also stored here.
 * \param v2 The divisor, it is left unchanged.
 * \return Error code.
 * 
 * Time complexity: O(n), the length of the vectors.
 */

int FUNCTION(igraph_vector,div)(TYPE(igraph_vector) *v1, 
				const TYPE(igraph_vector) *v2) {
  
  long int n1=FUNCTION(igraph_vector,size)(v1);
  long int n2=FUNCTION(igraph_vector,size)(v2);
  long int i;
  if (n1 != n2) {
    IGRAPH_ERROR("Vectors must have the same number of elements for swapping",
		 IGRAPH_EINVAL);
  }
  
  for (i=0; i<n1; i++) {
#ifdef DIV
    DIV(VECTOR(*v1)[i], VECTOR(*v1)[i], VECTOR(*v2)[i]);
#else
    VECTOR(*v1)[i] /= VECTOR(*v2)[i];
#endif
  }
  
  return 0;
}

#ifndef NOABS

int FUNCTION(igraph_vector,abs)(TYPE(igraph_vector) *v) {
#ifdef UNSIGNED
  /* Nothing do to, unsigned type */
#else
  long int i, n=FUNCTION(igraph_vector,size)(v);
  for (i=0; i<n; i++) {
    VECTOR(*v)[i]= VECTOR(*v)[i] >= 0 ? VECTOR(*v)[i] : -VECTOR(*v)[i];
  }
#endif
  
  return 0;
}

#endif

#ifndef NOTORDERED

/**
 * \function igraph_vector_minmax
 * \brief Minimum and maximum elements of a vector.
 * 
 * Handy if you want to have both the smallest and largest element of
 * a vector. The vector is only traversed once. The vector must by non-empty.
 * \param v The input vector. It must contain at least one element.
 * \param min Pointer to a base type variable, the minimum is stored
 *     here.
 * \param max Pointer to a base type variable, the maximum is stored
 *     here.
 * \return Error code.
 * 
 * Time complexity: O(n), the number of elements.
 */

int FUNCTION(igraph_vector,minmax)(const TYPE(igraph_vector) *v,
				   BASE *min, BASE *max) {
  long int n=FUNCTION(igraph_vector,size)(v);
  long int i;
  *min=*max=VECTOR(*v)[0];
  for (i=1; i<n; i++) {
    BASE tmp=VECTOR(*v)[i];
    if (tmp > *max) { 
      *max=tmp;
    } else if (tmp < *min) {
      *min=tmp;
    }
  }
  return 0;
}

/**
 * \function igraph_vector_which_minmax
 * \brief Index of the minimum and maximum elements
 * 
 * Handy if you need the indices of the smallest and largest
 * elements. The vector is traversed only once. The vector must to
 * non-empty.
 * \param v The input vector. It must contain at least one element.
 * \param which_min The index of the minimum element will be stored
 *   here.
 * \param which_max The index of the maximum element will be stored
 *   here.
 * \return Error code.
 * 
 * Time complexity: O(n), the number of elements.
 */

int FUNCTION(igraph_vector,which_minmax)(const TYPE(igraph_vector) *v,
					 long int *which_min, long int *which_max) {

  long int n=FUNCTION(igraph_vector,size)(v);
  long int i;
  BASE min, max;
  *which_min=*which_max=0;
  min=max=VECTOR(*v)[0];
  for (i=1; i<n; i++) {
    BASE tmp=VECTOR(*v)[i];
    if (tmp > max) {
      max=tmp;
      *which_max=i;
    } else if (tmp < min) {
      min=tmp;
      *which_min=i;
    }
  }
  return 0;
}

#endif

/**
 * \function igraph_vector_isnull
 * \brief Are all elements zero?
 * 
 * Checks whether all elements of a vector are zero.
 * \param v The input vector
 * \return Boolean, \c TRUE if the vector contains only zeros, \c
 *    FALSE otherwise.
 * 
 * Time complexity: O(n), the number of elements.
 */

igraph_bool_t FUNCTION(igraph_vector,isnull)(const TYPE(igraph_vector) *v) {

  long int n=FUNCTION(igraph_vector,size)(v);
  long int i=0;
  
#ifdef EQ
  while (i<n && EQ(VECTOR(*v)[i], ZERO)) {
#else
  while (i<n && VECTOR(*v)[i]==ZERO) {
#endif
    i++;
  }
  
  return i==n;
}

#ifndef NOTORDERED

int FUNCTION(igraph_i_vector,intersect_sorted)(
    const TYPE(igraph_vector) *v1, long int begin1, long int end1,
    const TYPE(igraph_vector) *v2, long int begin2, long int end2,
    TYPE(igraph_vector) *result);

/**
 * \function igraph_vector_intersect_sorted
 * \brief Calculates the intersection of two sorted vectors
 *
 * The elements that are contained in both vectors are stored in the result
 * vector. All three vectors must be initialized.
 *
 * </para><para>
 * Instead of the naive intersection which takes O(n), this function uses
 * the set intersection method of Ricardo Baeza-Yates, which is more efficient
 * when one of the vectors is significantly smaller than the other, and
 * gives similar performance on average when the two vectors are equal.
 *
 * </para><para>
 * The algorithm keeps the multiplicities of the elements: if an element appears
 * k1 times in the first vector and k2 times in the second, the result
 * will include that element min(k1, k2) times.
 *
 * </para><para>
 * Reference: Baeza-Yates R: A fast set intersection algorithm for sorted
 * sequences. In: Lecture Notes in Computer Science, vol. 3109/2004, pp.
 * 400--408, 2004. Springer Berlin/Heidelberg. ISBN: 978-3-540-22341-2.
 *
 * \param v1 the first vector
 * \param v2 the second vector
 * \param result the result vector, which will also be sorted.
 *
 * Time complexity: O(m log(n)) where m is the size of the smaller vector
 * and n is the size of the larger one.
 */
int FUNCTION(igraph_vector,intersect_sorted)(const TYPE(igraph_vector) *v1,
    const TYPE(igraph_vector) *v2, TYPE(igraph_vector) *result) {
  long int size1, size2;

  size1 = FUNCTION(igraph_vector,size)(v1);
  size2 = FUNCTION(igraph_vector,size)(v2);

  FUNCTION(igraph_vector,clear)(result);

  if (size1 == 0 || size2 == 0)
    return 0;

  IGRAPH_CHECK(FUNCTION(igraph_i_vector,intersect_sorted)(
    v1, 0, size1, v2, 0, size2, result));
  return 0;
}

int FUNCTION(igraph_i_vector,intersect_sorted)(
    const TYPE(igraph_vector) *v1, long int begin1, long int end1,
    const TYPE(igraph_vector) *v2, long int begin2, long int end2,
    TYPE(igraph_vector) *result) {
  long int size1, size2, probe1, probe2;

  if (begin1 == end1 || begin2 == end2)
    return 0;

  size1 = end1 - begin1;
  size2 = end2 - begin2;

  if (size1 < size2) {
    probe1 = begin1 + (size1 >> 1);      /* pick the median element */
    FUNCTION(igraph_i_vector,binsearch_slice)(v2, VECTOR(*v1)[probe1], &probe2, begin2, end2);
    IGRAPH_CHECK(FUNCTION(igraph_i_vector,intersect_sorted)(
        v1, begin1, probe1, v2, begin2, probe2, result
    ));
    if (!(probe2 == end2 || VECTOR(*v1)[probe1] < VECTOR(*v2)[probe2])) {
      IGRAPH_CHECK(FUNCTION(igraph_vector,push_back)(result, VECTOR(*v2)[probe2]));
      probe2++;
    }
    IGRAPH_CHECK(FUNCTION(igraph_i_vector,intersect_sorted)(
        v1, probe1+1, end1, v2, probe2, end2, result
    ));
  } else {
    probe2 = begin2 + (size2 >> 1);      /* pick the median element */
    FUNCTION(igraph_i_vector,binsearch_slice)(v1, VECTOR(*v2)[probe2], &probe1, begin1, end1);
    IGRAPH_CHECK(FUNCTION(igraph_i_vector,intersect_sorted)(
        v1, begin1, probe1, v2, begin2, probe2, result
    ));
    if (!(probe1 == end1 || VECTOR(*v2)[probe2] < VECTOR(*v1)[probe1])) {
      IGRAPH_CHECK(FUNCTION(igraph_vector,push_back)(result, VECTOR(*v2)[probe2]));
      probe1++;
    }
    IGRAPH_CHECK(FUNCTION(igraph_i_vector,intersect_sorted)(
        v1, probe1, end1, v2, probe2+1, end2, result
    ));
  }

  return 0;
}

/**
 * \function igraph_vector_difference_sorted
 * \brief Calculates the difference between two sorted vectors (considered as sets)
 *
 * The elements that are contained in only the first vector but not the second are
 * stored in the result vector. All three vectors must be initialized.
 *
 * \param v1 the first vector
 * \param v2 the second vector
 * \param result the result vector
 */
int FUNCTION(igraph_vector,difference_sorted)(const TYPE(igraph_vector) *v1,
    const TYPE(igraph_vector) *v2, TYPE(igraph_vector) *result) {
  long int i, j, i0, j0;
  i0 = FUNCTION(igraph_vector,size)(v1);
  j0 = FUNCTION(igraph_vector,size)(v2);
  i = j = 0;

  if (i0 == 0) {
    /* v1 is empty, this is easy */
    FUNCTION(igraph_vector,clear)(result);
    return IGRAPH_SUCCESS;
  }

  if (j0 == 0) {
    /* v2 is empty, this is easy */
    IGRAPH_CHECK(FUNCTION(igraph_vector,resize)(result, i0));
    memcpy(result->stor_begin, v1->stor_begin, sizeof(BASE) * (size_t) i0);
    return IGRAPH_SUCCESS;
  }

  FUNCTION(igraph_vector,clear)(result);

  /* Copy the part of v1 that is less than the first element of v2 */
  while (i < i0 && VECTOR(*v1)[i] < VECTOR(*v2)[j]) { i++; }
  if (i > 0) {
    IGRAPH_CHECK(FUNCTION(igraph_vector,resize)(result, i));
    memcpy(result->stor_begin, v1->stor_begin, sizeof(BASE) * (size_t) i);
  }

  while (i < i0 && j < j0) {
    BASE element = VECTOR(*v1)[i];
    if (element == VECTOR(*v2)[j]) {
      i++; j++;
      while (i < i0 && VECTOR(*v1)[i] == element) { i++; }
      while (j < j0 && VECTOR(*v2)[j] == element) { j++; }
    } else if (element < VECTOR(*v2)[j]) {
      IGRAPH_CHECK(FUNCTION(igraph_vector,push_back)(result, element));
      i++;
    }
    else j++;
  }
  if (i < i0) {
    long int oldsize = FUNCTION(igraph_vector,size)(result);
    IGRAPH_CHECK(FUNCTION(igraph_vector,resize)(result, oldsize+i0-i));
    memcpy(result->stor_begin+oldsize, v1->stor_begin+i, 
	   sizeof(BASE)* (size_t) (i0-i));
  }

  return 0;
}

#endif

#if defined(OUT_FORMAT)

#ifndef USING_R
int FUNCTION(igraph_vector,print)(const TYPE(igraph_vector) *v) {
  long int i, n=FUNCTION(igraph_vector,size)(v);
  if (n!=0) {
#ifdef PRINTFUNC
    PRINTFUNC(VECTOR(*v)[0]);
#else
    printf(OUT_FORMAT, VECTOR(*v)[0]);
#endif
  }
  for (i=1; i<n; i++) {
#ifdef PRINTFUNC
    putchar(' '); PRINTFUNC(VECTOR(*v)[i]);
#else
    printf(" " OUT_FORMAT, VECTOR(*v)[i]);
#endif
  }
  printf("\n");
  return 0;
}

int FUNCTION(igraph_vector,printf)(const TYPE(igraph_vector) *v,
																	 const char *format) {
  long int i, n=FUNCTION(igraph_vector,size)(v);
  if (n!=0) {
    printf(format, VECTOR(*v)[0]);
  }
  for (i=1; i<n; i++) {
    putchar(' '); printf(format, VECTOR(*v)[i]);
  }
  printf("\n");
  return 0;
}

#endif

int FUNCTION(igraph_vector,fprint)(const TYPE(igraph_vector) *v, FILE *file) {
  long int i, n=FUNCTION(igraph_vector,size)(v);
  if (n!=0) {
#ifdef FPRINTFUNC
    FPRINTFUNC(file, VECTOR(*v)[0]);
#else
    fprintf(file, OUT_FORMAT, VECTOR(*v)[0]);
#endif
  }
  for (i=1; i<n; i++) {
#ifdef FPRINTFUNC
    fputc(' ', file); FPRINTFUNC(file, VECTOR(*v)[i]);
#else
    fprintf(file, " " OUT_FORMAT, VECTOR(*v)[i]);
#endif
  }
  fprintf(file, "\n");
  return 0;
}

#endif

int FUNCTION(igraph_vector,index)(const TYPE(igraph_vector) *v,
                                  TYPE(igraph_vector) *newv,
                                  const igraph_vector_t *idx) {
  
  long int i, newlen=igraph_vector_size(idx);
  IGRAPH_CHECK(FUNCTION(igraph_vector,resize)(newv, newlen));
  
  for (i=0; i<newlen; i++) {
    long int j=(long int) VECTOR(*idx)[i];
    VECTOR(*newv)[i] = VECTOR(*v)[j];
  }
  
  return 0;
}
